function [gz]=grav1D(H,rho,Lx,Ly)

n_slabs=length(H); %number of slabs
nx=length(rho(1,:)); %number of points in x direction
ny=1; %points in y direction

%wave numbers
kx=0:1:nx-1;
ky=0:1:ny-1;
[kxi,kyi]=meshgrid(kx,ky);

G=6.673e-11; %gravity constant

%calculate Omega
O=(pi*kxi/Lx).^2+(pi*kyi/Ly).^2;
SqO=O.^0.5;

%in order to prevent division by zero:
SqO(1,1)=1;

gz=zeros(ny,nx);
z=0; %elevation above the top slap is zero (z=0 means that we stand on it)

for k=1:length(H), %loop over all slabs

    %fft of mass distribution
    R=fct2c(rho(k,:));
    
    %z is the depth to the top of slab no k;
   
    F=2*pi*G*(exp(SqO*H(k))-1).*exp(-SqO*z).*R./SqO;
    
    %wave no k=0 and l=0, stored in (1,1), require special treatment
    F(1,1)=2*pi*G*R(1,1)*H(k);
    
    %contribution from slab no k to gravity field
    gz=gz+ifct2c(F);
        
    %update depth
    z=z+H(k);
    
end;


